Rare regions of the SIS model on Barabasi- Albert networks 



O 
(N 



in 

"a 

C/2 



Geza Odor 

Research Centre for Natural Sciences, Hungarian Academy of Sciences, 
MTA TTK MFA, P. O. Box 49, H-1525 Budapest, Hungary 
(Dated: March 6, 2013) 

I extend a previous work to susceptible-infected-susceptible (SIS) models on weighted Barabasi- 
Albert scale-free trees. Numerical evidence is provided that phases with slow, power-law dynamics 
emerge as the consequence of quenched disorder in topologies previously studied with Contact 
Process. I compare simulation results with a spectral analysis of the networks and show that the 
quenched mean-field (QMF) approximation provides a reliable, relatively fast method to explore 
activity clustering. This suggests that QMF can be used for describing rare-region effects and 
smeared phase transitions due to network inhomogeneities. Finite size study of the QMF shows the 
expected disappearance of the epidemic threshold Ac in the thermodynamic limit and an inverse 
participation ratio ~ 0.25, meaning localization in case of assortative weights. Contrary, for the 
multiplicative weights and the unweighted trees this value vanishes in the thermodynamic limit, 
suggesting only weak rare region effects in agreement with the dynamical simulation results. Strong 
corrections to the mean-field behavior in case of assortative weights explains the concave shape of 
the order parameter p(A) at the transition point. Application of this method to other models may 
reveal interesting rare-region effects, Griffiths Phases as the consequence of quenched topological 
heterogenities. 
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I. INTRODUCTION 

The research of nonequilibrium models has been a cen- 
tral research topic of statistical mechanics [H-Q . A fun- 
damental dynamical system model to understand them 
is the Contact Process (CP) [1, in which sites can be 
either occupied (infected) or empty (susceptible). By de- 
creasing the infection rate of the neighbors A/fc, where k 
is the degree of the vertex, a continuous phase transition 
occurs at some Ac critical point from active to inactive 
steady state. The latter is also called absorbing, because 
no spontaneous activation of sites is allowed and the den- 
sity of infection (p) is zero. 

Recently the interest has been shifted from models, 
defined on Euclidean, regular lattices to processes living 
on general networks @, 0] • The effects of the heteroge- 
neous topological structures [1, Q is not yet fully under- 
stood. In particular in the case of ubiquitous scale-free 
(SF) networks exhibiting P{k) ~ fc-T de gree distri- 
bution of the nodes @, 0] the location of the phase tran- 
sition and the singular behavior around is still a debated 
issue. Numerical simulations and theoretical ap- 

proaches based on the heterogeneous mean-field (HMF) 
theory [Hllillili show strong effects of the network het- 
erogeneity on the behavior of the CP defined on complex 
networks. 

Although SF-s exhibit infinite topological dimension 
(d), defined by A*" cx r'', where N is the number of nodes 
within the (chemical) distance r, simple mean- field ap- 
proximations cannot capture several important features. 
Studies of the CP, as well as other processes [1, Q have 
shown that quenched disorder in networks is relevant in 
the dynamical systems defined on top of them. Very 
recently it has been shown, [l6l - [l8| that generic slow 



(power-law, or logarithmic) dynamics is observable by 
simulating CP on networks with finite d. This obser- 
vation is relevant for recent developments in dynamical 
processes on complex networks such as the sirnple model 
of "working memory" [loj . brain dynamics [20l. social 
networks with heterogeneous communities [2l| . or slow 
relaxation in glassy systems p^ . 

Slow dynamics can be the consequence of bursty be- 
havior of agents connected by small world networks re- 
sulting in memory effects |23| . An independent cause 
is related to such arbitrarily large {I < N), correlated 
rare-regions (RR), which have long lifetime r cx exp(Z) in 
the active phase (A > A°). In the re gion < A < Ac, 
a so called Griffiths Phase (CP) (23. l25| develops with 
an algebraic density decay p cx a being a non- 

universal exponent. This can only be understood by non- 
perturbative methods j26l - [29| . 

More recently the possibility of such slow dynamics has 
been investigated on BA networks with 7 = 3 possessing 
infinite d [30|. Very extensive simulations showed linear 
leading order density decay and p{X,t — >■ oo) cx |A — Ad 
behavior with logarithmic corrections in agreement with 
the HMF approximations. It was also pointed out that 
in case of the loop-less BA trees the epidemic propaga- 
tion slows down, and a nontrivial critical density decay 
emerges with p{t, Ac) cx a ~ 0.5. Furthermore, when 
k dependent weighting was also applied, suppressing hubs 
or making the network disassortatative GP-like regions 
could be observed in the simulations. A systematic finite 
scaling study revealed that these power-laws saturate in 
the N ^ oo thermodynamic limit, suggesting smeared 
phase transitions [2^. This can be understood because 
even infinite dimensional RR-s can be embedded in net- 
works with d = oo causing finite p(t — > oo,N — > oo). 
A numerical percolation analysis has strengthened this 
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view indeed. 

I continue this work [s^l to see if quenched mean- 
field (QMF) theory with spectral decomposition (SD) can 
help to understand rare-region effects in complex net- 
works. For this reason I study the susceptible- infected- 
susceptible (SIS) model [3l|, instead of the CP because it 
possesses real symmetric adjacency matrix with positive 
eigenvalues. The susccptible-infcctcd-susceptiblc (SIS) 
model is a two state system, in which infected sites prop- 
agate the epidemic (venerate all neighbors) with rate A, 
or recover with rate 1. 

This extension is far from trivial, because the epidemic 
transmitting capability of infected nodes is higher, i.e. 
it is not normalized by the number of outgoing edges. 
Therefore the emergence of localized rare-regions (RR) 
is questionable. I provide extensive simulation results 
showing numerical evidence for GP like regions with slow 
dynamics similarly as in case of the CP. The QMF and 
SD analysis set up for these cases is in good agreement 
with the dynamical simulations. 

II. NETWORK MODELS 



where the node numbers i and j correspond to the time 
step when they were connected the network. Since the 
degree of nodes decreases as ki oc {N/iY^'^ during this 
process, this selection with x > favors connection be- 
tween unlike nodes and suppresses interactions between 
similar ones. In [30l | simulations showed evidence for a 
phase with power-law dynamics of the CP for a; = 2, 3 
with CP. This paper is concerned about x = 2 networks. 

The presence of these weights affects the dynamics of 
the SIS. Thus, the rate at which a healthy vertex i be- 
comes ill on contact with an infected (active) vertex j is 
proportional to Awy , therefore the epidemic can in princi- 
ple become trapped in isolated connected subsets. These 
regions of size I are rare in general: P{1) oc cxp(— Z), but 
can exhibit exponentially long lifetimes t{1) cx cxp(Z). 
In the healthy (inactive) phase they provide the leading 
order contribution to the density of infected sites 

p{t) - J IP{1) exp{-t/T)dl , (3) 

which in the saddle point approximation results in p{t) ^ 
t-°' decay 



I consider here SIS on BA networks in particular for 
loop-less and weighted cases as described in [sO]- This 
permits very simple and fast construction, in contrast 
with other standard network generation models, e.g fs^ . 
The BA growth starts with a fully connected graph of 
size No = 10 nodes, but comparisons with Nq = 5 and 
20 have also been done to test any dependence. For BA 
at each time step s, a new vertex with m edges is added to 
the network and connected to an existing vertex s' of de- 
gree ks' with probability IIs-j.^' = k^i / ^^ksn . This 
process is iterated until reaching the desired network size 
N . The resulting network has a SF degree distribution 
P{k) ~ For m = 1 we obtain a BA tree (BAT) 

topology, while for the looped case I applied to = 3. 

Binary (non-weighted) BA networks can be trans- 
formed into weighted ones by assigning to every edge 
connecting vertexes i and j a symmetric weight uJij. In 
(30j two different network topology dependent weight as- 
signment strategy was introduced in order to slow down 
and localize epidemics. 

(i) Weighted BA tree I (WBAT-I): Muhiphcativc 
weights, suppressing the infection capability of highly 
connected nodes 



(1) 



where wq is an arbitrary scale and u is a characteristic 
exponent with v > 0. This can model internal limitations 
of hubs, like the sub-linear Heap's law [s^. In this paper 
I study the case with ^/ = 1.5 exponent. 

(ii) Weighted BA tree II (WBAT-II): Disassortatative 
weighting scheme according to the age of nodes in the 
network construction 



N 



(2) 



III. SPECTRAL ANALYSIS 

In [s^l the heterogeneous mean-field (HMF) analysis of 
CP was worked out for these network models. However, 
extensive simulations showed different dynamical behav- 
iors, except from the looped BA CP model case. Since 
HMF can't take into account rare region effects, nor mod- 
els on trees that conclusion was not very surprising. 

A mean-field theory, capable of describing the effects 
of quenched topologies of the network on which SIS is de- 
fined is expected to give better agreement with numerical 
simulations. The quenched mean-field (QMF) approach 
is based on the rate equation for Pi{t), the infection prob- 
ability of node i at time t (33| : 



dpijt) 
dt 



^-p,it) + il-p.,{t))Y,A^,XpJit) , (4) 



where Aij is an element of the adjacency matrix assigned 
with 1, if there is an edge between nodes i and j or 
otherwise. This equation can be generalized by replacing 
the adjacency matrix with the weighted adjacency matrix 
Bij = AijCUij, having weighted elements Bij G [0, 1] due 
to its real symmetry. 

For t ^ oo the system evolves into a steady state, with 
the probabilities expressed as 



1 



^lljBijPj 



(5) 



Stability analysis shows that pi > above a Ac epidemic 
threshold, with finite order parameter (prevalence) p = 

{P^) ■ 
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In the SD approach one expands pi in the space of 
eigenvectors of the adjacency matrix Aij (or Bij for 
weighted case) as 



(6) 



This extension can be done for real and symmetric 
weights, when the eigenvectors /(A) span a complete 
orthonormal basis. The Perron-Frobcnius theorem as- 
serts that a real square matrix with positive entries has a 
unique largest real eigenvalue Kmax = Ai > A2 > . . . h-N 
and that the corresponding eigenvector /(Ai) has strictly 
positive components. In this basis Eq. ([5]) can be ex- 
pressed by the coefficients c(A) 



N 



:(A) = aEA'c(A')E 



/,(A)/,(A') 



tl + AEAAc(A)/,(A) 



(7) 



This gives Ac = 1/Ai for the epidemic threshold and in its 
neighborhood the order parameter can be approximated 
by the eigenvectors of the largest eigenvalues 



p{X) ssaiA + a2A2 + ... , 
where A = AAi — 1^1 with the coefficients 



N 



N 



E/,(A,)/[ivE/f(A,)]. 



(8) 



(9) 



i=l 



i=l 



A homogeneous solution implies that a finite fraction of 
vcrtcxcs arc infected right above Ac and ai is the or- 
der of 0(1). That would mean that the components of 
/(Ai) are localized. On the other hand quenched inho- 
mogeneous topologies can also imply inhomogeneous pi 
distributions and one can assume that they result in rare 
region effects, as in 

To describe localization in the components of /(Ai) 
[33} suggested to use the inverse participation ratio 



N 



/Pi?(A)^E//(A), 



(10) 



which in the limit TV — >■ cx) is of the order of 0(1) when 
the eigenvector is localized. Contrary, when IPR{A) — > 
this state is delocalized. Eq. © implies that for localized 
principal eigenvector ai ~ 0(1/A^), thus p « aiA ^ 
0{1/N), the disease is localized on a finite number Np 
of vertexes. On the other hand, when /(Ai) is delocalized 
the disease infects a finite fraction of vertexes for A > Ac- 
Goltscv et al. [s^ analyzed artificial and real SF net- 
works and found that in case of localized cases the epi- 
demic threshold was actually absent and a real epidemic 
affecting a finite fraction of vertexes occurred after a 
smooth crossover at higher values of A. This is in agree- 
ment with the smeared phase transition scenario pro- 
posed to explain the numerical results for CP on weighed 
BA trees in [sO]. In those systems rare- regions effects 



seemed to arise, causing power-law density decays, which 
ultimately crossed over to saturation to finite p-s in the 
N ^ 00 limit. Right above Ac similar concave shaped 
p{X) was obtained as in the localized cases of [s^. Here 
I investigate the situation for SIS instead of CP models, 
to see if this relation holds on. This has been done using 
finite size scaling analysis of the quantities: A = 1/Ai, 
IPR and ai for i = 1, 2, 3. 

It has been proven [34| that for random, unweighted 
SF networks, with power-law degree distribution P{k) oc 
k^'^ the maximal eigenvalue scales as Ai oc V kmax for 
7 > 2.5. Furthermore, the maximal degree of the net- 
work satisfies k^ax = min [TV^/^, Afi/(7-i)] ^ due to the 
structural cutoff of a finite network with 7 < 3. There- 
fore, Al should scale with the network size as 



Ai(7V) cx TV^/" 



(11) 



for 7 = 3 considered here. 

Using the software package OCTAVE I generated the 
Bij matrices of BA, BAT, WBAT-I and WBAT-II net- 
works for several sizes up to iV = 200.000 and calculated 
the three largest eigenvalues and the corresponding eigen- 
vectors. From these I deduced 1/Ai, IPR and ai. The 
whole SD analysis was done using the sparse matrix func- 
tions of OCTAVE to handle Bij-s of the networks within 
reasonable computing times. For the largest sizes 1-2 
weeks of a CPU time was needed. Least-squares fitting 
with the form 



1/Ai = a + &(l/iV)^ 



(12) 



has been applied for the largest eigenvalues. As Table. |T] 
shows a good agreement was obtained with the finite size 
scaling expectation lim^v^oo = 1/Ai = of the QMF 
method. The N ^ 00 extrapolated critical threshold 
values converge to zero using the three parameter fit- 
ting form (|12p . The fitted power c agrees with exponent 
of (fTT|) reasonably well, except from the WBAT-I case, 
where QMF results for only < 6000 could be achieved. 

The IPR values decrease with 1/A^ and remain very 
small for the unweighted BA {IPR < 0.02 ) and BAT 
{IPR < 0.06) networks, albeit in the latter case the 
tendency seems to change for N > 32000 as shown on 
Fig. [TJ This suggests an agreement with the anomalous 
scaling behavior of CP on the BA tree reported in [30j . 
Here some clustering may be expected, which should be 
checked further by studying networks with N > 10^ 
nodes. Unfortunately this is out of the scope of the 
present study using OCTAVE, but is in preparation using 
graphics cards. For the WBAT-II case the IPR stabilizes 
to ~ 0.25(1) (see Fig. [5]), suggesting a strong epidemic 
localization. 

One can think that the initial compact seeds, from 
which the BA graphs originate can infiuence the clus- 
tering behavior of the epidemic in the steady state. By 
varying the seed size: No = 5, 10, 20 no measurable effect 
was found on the WBAT results. For the unweighted BA 
and BAT networks on can see differences between the fi- 
nite size results of the A'o = 20 and Nq = 10 cases (see 
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FIG. 1: (Color online) Finite size scaling of QMF SD results 
for BAT model for iV = 4000, 8000, 16000, 32000, 64000, 
100000. Filled symbols correspond to A*'o = 10, hollow ones 
to A^o ~ 20. Dashed line: least-squares fitting with the form 



TABLE I: Spectral QMF analysis resuhs for SIS in different 
networks 



Network 


1/Ai 


c 


IPR 


Ql 


ce2 


Q3 


BA 


0.017 


0.32 


0.014 


4 X 10""* 




10"* 


BAT 


~ 


0.24 


0.055 


10"^ 


10-' 


io-» 


WBAT-I 


~ 


0.9(1) 


~ 


0.35(1) 


5 X 10"* 


10-*^ 


WBAT-II 


0.001 


0.227 


0.25 


4 X 10"^ 


3 X lO"'^ 


10"^ 



Fig. [T]), but the iV — > oo asymptotic behavior appears to 
be the same. 

For the unweighted BA and BAT networks the eigen- 
value analysis results in ai » ai » for all sizes 
in agreement with j3 = 1 expectations for p slightly 
above the critical threshold. In the oo limit 

ai(iV) ~ (l/A^)"-^ and a2 — as ~ already for small 
sizes. For the weighted WBAT-I case ai remains con- 
stant ~ 0.35(1) while a2 ^ ^ (l/N). In case of the 
WBAT-II tree each ai is the same order of magnitude 
and vanish linearly with meaning a stronger cor- 

rection to the leading order linear scaling as one can see 
in Fig.[2j When one plots p(A) with these values one gets 
a concave curve from above and a tangential approach to 
Ac- Such steady state behavior has already been seen by 
simulations of CP on weighted B A networks [30| . In the 
next section I compare these QMF results with simula- 
tions of SIS on WBAT-I and WBAT-II networks. 



IV. SIS MODEL SIMULATIONS ON 
WEIGHTED TREES 

In simulations I considered the SIS model in contin- 
uous time as in (36j . to be in accordance with the rate 



FIG. 2: (Color online) Finite size scaling of QMF SD results 
for WBAT-II model for = 2000, 4000, ... 200000. Solid 
line: least-squares fitting for with the form (|12p . 

equations. At each time step a randomly chosen infected 
node recovers with probability ni/{ni + Ann), where rii 
is the the number of infected nodes and n„ is the to- 
tal number of links emanating from them. Complemen- 
tary, one of its randomly selected neighbor is infected, 
with probability Xni/{ni + An„). Following the reaction 
Ni and Nn are updated and the time is incremented by 
At = l/{ni + Xnn)- The time is measured by these Monte 
Carlo steps (MCs) and shown to be dimensionless on the 
figures. These processes are iterated until t < t„iax, or 
until the epidemic stops {Ni = 0). 

The networks were generated via the BA linear pref- 
erential attachment rule following an initial fully 
connected seed of A''o = 20. Neighbor indices of sites 
are stored in a dense matrix to save memory, thus up 
to iV = 6 X 10^ sized networks could be studied. The 
initial state was fully active and the concentration of in- 
fected sites was followed up to t^ax = 4 x 10^ MCs. 
Density decay runs were repeated and averaged over 
~ lO** and up to 10^ independent network realizations 
for WBAT-I and WBAT-II, respectively. The behavior 
in the active steady state was also investigated by mea- 
suring p{X,t — >■ oo) deeply in the saturation region. 

Figure [SI suggests that p{t) of the WBAT-I model ex- 
hibits A dependent power-laws for t > 3 x 10^ MCs. This 
can be observed in networks with N = 2 x 10^ nodes 
in the region 37.72 < A < 37.77. The final slope of at 
the lowest, A = 37.72 power-law curve is: a ~ 0.343. 
By increasing the system size a level-off of these curves 
can be observed as in [30| . However, this range is very 
narrow both in A and the effective exponent a. Further- 
more, curvature can also be seen on the log. -log. plots 
on some of them in this 'scaling' region, which makes the 
power-law dynamics assumption questionable. In com- 
parison, for CP in the region 140 < A < 145 power-laws 
with 0.3 < a < 1 can be seen clearly [3^. Therefore, 
rare-region effects seem to be much weaker in agreement 
with QMF SD analysis. When we plot the same data on 
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FIG. 3: (Color online) Density decay as a function of time 
for the SIS model on weighted BA trees generated with the 
WBAT-I scheme with exponent i/ = 1.5. Network size A'^ = 
2x10'^. Different curves correspond to A = 37.72, 37.74, 37.75, 
37.76, 37.77 (from bottom to top). Inset: The same data 
plotted on the ln{— \n{p{t))) vs. \n{t) scale (with parameters 
from top to bottom order). 
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FIG. 4: (Color online) Density decay as a function of time 
for the SIS on weighted BA trees generated with the age- 
dependent (disassortative) WBAT-II scheme with exponent 
X = 2. Network size N = 10^. Different curves correspond to 
A = 1.657, 1.66, 1.661, 1.663, 1.67, 1.68, 1.69 (from bottom 
to top). Dashed line; power-law fit. Inset; Steady state 
density (bullets) above the epidemic threshold. The curve 
shows power-law fitting. 



the ln(— ln(p(t))) vs. ln(i) scale, straight lines appear, 
suggesting stretched exponential behavior (see inset of 
Figure El). 

In the case of WBAT-II network the power-laws decays 
of p{t, A) appears much earlier, already or i > 50.000 
MCs and more pronouncedly, albeit corrections to scaling 
are non-negligible again. The final slope of the lowest, 
A = 1.657 curve is: a ~ 0.674, obtained by fitting for 
t > 20.000 MCs, which is far away from the HMF critical 



exponent: a = 1. This behavior is located in the 1.657 < 
A < 1.67 region, by much lower values than those of the 
CP dOl- A smeared transition is expected again, since 
the densities in the steady state p{t — >■ oo) increase with 
N. As the inset of Fig.|4]shows p(A, t — )- oo) can be fitted 
using 

p{X,t ^ oo) =C{X- Xcf (13) 

with an order parameter exponent /3 ~ 1.22(1) near Ac — 
1.64(1) for networks with size iV 4 x 10*^. 

V. DISCUSSION AND CONCLUSIONS 

Understanding effects of heterogenities in nonequilib- 
rium, dynamical systems is a challenging open field. 
References [l^ [l^ concluded that finite topological di- 
mension is a necessary condition for observing Griffiths 
Phases and activated scaling in case of the basic model of 
nonequilibrium system, the Contact Process. In case of 
CP on certain weighed networks numerical evidence was 
shown for generic power-laws, but in the thermodynamic 
limit this phase seems to disappear and a smeared phase 
transition exists, due to the infinite dimensional corre- 
lated rare- regions [13] ■ In contrast with these strong dis- 
order renormalization studies of the random transverse- 
field Ising model found Griffiths singularities (ioj even in 
the infinite dimensional Erdos-Renyi random graphs [4l| . 

Very recently Griffiths phases have been reported in 
a study of the Random Transverse Ising Model on com- 
plex networks with a scale-free degree distribution regu- 
larized by and exponential mtoS p{k) oc k^'^ exp[—k/S,] 
[38| . This model was devised to understand the relation 
between the onset of the superconducting state with the 
particular optimum heterogeneity in granular supercon- 
ductors. On quenched networks a phase transition at 
zero temperature and a Griffith phase with decreasing 
size as the function of the cutoff ^ was found. The sce- 
nario is similar to our case, because increasing the value 
of the cutoff is like going to the thermodynamic limit of 
the SF network. Another fresh model study [s^ argues 
again for the existence of CP in the SIS model defined 
on unclustered, deterministic SF networks, with 7 > 3, 
below the percolation threshold. 

In this work I provided spectral decomposition and 
QMF approximation to SIS models on different scale-free 
networks. This analysis has been supplemented with ex- 
tensive numerical simulations showing dynamical effects 
of the topological disorder. Activity localization of the 
eigenvectors characterized by the IPR number predicts 
strong rare-region effects in the case of weighted WBAT- 
II trees. Numerical simulations exhibit a A parameter 
region, with continuously changing power-laws decays. 
Still, in the N ^ 00 limit the p{t) curves saturate, as 
for smeared phase transitions, discussed in (30j . For the 
WBAT-I tree the picture is less clear. There is a narrow 
CP like region at very late times {t > 3 x 10^ MCs), but 
the p{t) curves exhibit strong corrections and a stretched 
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exponential dynamics looks more reasonable. This see- 
nario is strengthened by the QMF SD analysis, which 
results in a vanishing IPR, suggesting only weak rare- 
region effects here. 

This study goes one step further than [s^ by analyzing 
spectral data of different graphs via finite size scaling. In 
particular numerical evidence is shown for the disappear- 
ance of Ac with the expected scaling law in the iV — > oo 
limit. Fluctuations, omitted by the QMF approximation 
can be relevant, indeed Ac > obtained by simulations 
of SIS on finite graphs. However, finite size study sug- 
gest Ac = in all cases, in accordance with smeared 
phase transition. Still the QMF SD analysis seems to 
be a relatively fast, promising way to explore topological 
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